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Abstract 

In recent years, a number of experimental studies have been conducted to investi- 
gate the mechanical behavior and damage mechanisms of articular cartilage under 
impact loading. Some experimentally observed results have been explained using a 
non-linear viscoelastic impact model. At the same time, there is the need of simple 
mathematical models, which allow comparing experimental results obtained in drop 
impact testing with impact loads of different weights and incident velocities. The 
objective of this study was to investigate theoretically whether the main features of 
articular impact could be qualitatively predicted using a linear viscoelastic theory 
or the linear biphasic theory. In the present paper, exact analytical solutions are 
obtained for the main parameters of the Kelvin- Voigt and Maxwell impact models. 
Perturbation analysis of the impact process according to the standard viscoelastic 
solid model is performed. Asymptotic solutions are obtained for the drop weight 
impact test. The dependence of the coefficient of restitution on the impact or pa- 
rameters has been studied in detail. 
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Nomenclature 
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damper coefficient 


D 


discriminant of the characteristic equation 
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incremental dynamic modulus 


E 


maximum incremental dynamic modulus 




modulus at stresses of 10 MPa 




coefficient of restitution 


F 


contact force 


Fm 


maximum contact force 


9 


gravitational acceleration 


h 


cartilage layer thickness 


h 


drop height of the impactor 




aggregate modulus 


k 


stiffness coefficient 


k u k 2 


spring stiffnesses in the standard solid model 


k 


instantaneous stiffness 


h 


long-term stiffness 


m 


impactor mass 


t 


time variable 


to 


impact duration 




time to maximum displacement 






vo 


initial impact velocity 


X 


displacement 


X 


velocity 


X 


acceleration 




maximum displacement 



f3 damping coefficient in the Kelvin- Voigt model 

Pi real part of complex roots of the characteristic equation 

Am percentage increase in mass of cartilage sample 

e strain 

So non-dimensional parameter accounting for the gravitational effect 

( loss factor in the Maxwell model 

(i imaginary part of complex roots of the characteristic equation 

rj loss factor in the Kelvin- Voigt model 

k cartilage permeability 

xi, x 2 spring stiffnesses in the standard solid model 

A Lame coefficient 

Ai root of the characteristic equation 

A non-dimensional parameter in the standard solid model 

fi Lame coefficient 

£ non-dimensional displacement 

p ratio of the long-term and instantaneous stiffnesses 

a stress 

r non-dimensional time 

td typical diffusion time 

tr relaxation time 

\?(t) dimensionless relaxation function 

uj angular frequency of damped oscillations 

ujo angular frequency of undamped oscillations 

1 Introduction 

Articular cartilage is a soft hydrated tissue covering the end of each bone at the joints. 
Cartilage has no known function other than maintaining mechanical competence of joints, 
allowing bones to move against one another without friction. But there is no need to 
underline its significance to health of a human body, since almost all the load transmitted 
by a human joint goes through the articular cartilage, and it prevents biomechanical 
damage caused by severe loading including impact loading. It is believed that severe 
articular impact can initiate post-traumatic arthritis [T1I2] . An impact loading of the joint 
constitutes the action of extremely high non-physiological loads applied very rapidly (for 



instance, due to a car accident, sports injury, or a fall from a height). 



In recent years, a number of experimental studies have been conducted to investigate the 
mechanical behavior and damage mechanisms of articular cartilage under impact loading 
J. In particular, the experimental data on relative dissipation of the impact energy 



AE/Eq versus overall impactor energy E obtained in [6] were fitted with quadratic 
curves. Here, E = AE = m(v\ — v$)/2, v and Vi are the initial impact and 

rebound velocities, respectively, m is the impactor mass. Since, v\ = — e*^o, where e* is the 
coefficient of restitution, we easily get AE/E = 1 — e*. Thus, the experimental data and 
fitting curves for dissipation of the impact energy [6] can be recalculated in terms of the 
coefficient of restitution as presented in Fig. [TJ which shows a non-monotonic dependence 
of e* on Vq. Some experimentally observed results have been explained using a non-linear 
viscoelastic impact model [7j. At the same time, there is the need of a simple mathematical 
model, which allows comparing experimental results obtained in drop impact testing with 
impact loads of different weights and incident velocities. 
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Fig. 1. Coefficient of restitution e* versus the impact velocity vo for articular cartilage samples 
of different thicknesses. Based on the experimental data and fitting curves obtained in [6]. 

A variety of mathematical models were suggested to describe the stress-strain response of 
articular cartilage that represents a multiphasic, structurally complex material possess- 
ing viscoelastic properties. It is long known that articular cartilage possesses viscoelastic 
properties [8ll9] , though there is no direct correspondence between viscoelastic parameters 
and parameters of the biphasic/poroelastic models of cartilage. The biphasic theory [ID] , 
which models the tissue as a mixture of a solid phase and a fluid phase, has demonstrated 
very good agreement with experimental results in the creep and stress relaxation tests 
[IT] . The objective of this study was to investigate theoretically whether the main fea- 
tures of articular impact observed in [6117] could be qualitatively predicted using a linear 
viscoelastic theory or the linear biphasic theory. 

The rest of the paper is organized as follows. In Sections [2] and [3| we consider in detail the 
viscoelastic Kelvin- Voigt and Maxwell impact models, respectively. Since some elements 
of the presented solutions are known in the literature, we pay a particular attention to 
the evaluation of the contact force, and impactor displacement, at the time 

moments £m and i m , when the force and displacement reach their maxima, Fm and x m) 
respectively. In Section [4j we outline a closed form solution of the impact equation in the 



case of standard solid model. In order to get analytical approximations, we consider the 
standard solid model as a perturbation of the Kelvin- Voigt (Section [5]) or the Maxwell 
model (Section [6]). In particular, simple analytical approximations are derived for the 
impact duration, i c , and for the coefficient of restitution, e*. In Sections [7] and [8| we 
consider the influence of the gravity effect on these parameters in the framework of the 
Kelvin- Voigt and Maxwell models for drop weight impact. In Section [9j we develop an 
asymptotic model for the force-displacement relationship in the indentation problem for a 
thin biphasic layer corresponding to the conditions of the so-called blunt impact, when the 
specimen thickness is much smaller than the radius of a flat-ended cylindrical impactor. An 
example of application of the developed linear theory of viscoelastic impact for analyzing 



experimental data is given in Section [TOj Finally, in Sections Y\_ and 12, we outline a 
discussion of the results obtained and formulate our conclusions. 



2 Viscoelastic Kelvin— Voigt impact model 



In this section, the deformation of articular cartilage layer is modeled schematically as a 
parallel combination of linear spring k and dashpot b (Fig. [2]). Dynamic balance between 
the force of cartilage reaction 

F = kx + bx (1) 

and the force of body inertia mx governs the development of collision. According to 
Newton's second law, the differential equation of the impact has the form 



mx + bx + kx = 0, t G [0, t c ] , 



(2) 



where t c is the contact duration, that is t c denotes the instant, when the cartilage reaction 
force changes its sign, or the indenter acceleration vanishes. 
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Fig. 2. Impact viscoelastic Kelvin- Voigt model. 
The initial conditions for Eq. ((2j) are as follows: 

x(0) = 0, x(0) = v . 



(3) 



The impact duration is determined by the condition 



kx + bx 



t=t c 



(4) 



or, in view of Eq. ([2]), by the condition 



x 



t=t c 



0. 



(5) 



The impact problem ([2]), ^ has the following well-known solution [T2] : 



x(t) = —e~ pt smut, t e [0,t e ]. 

UJ 



Here we used the notation 



UJ n 



m 



oj 2 = uj 2 -f3 2 , f3 = 



2m 



(6) 



(7) 



We assume that uj > f3. 



Fig. [3] shows the behavior of the dimensionless quantities ujqx/vq, F/(mvQOJo), and x/vq 
with respect to time. Observe that the time moment £m, when the contact force reaches 
its maximum, approaches the initial moment of impact as the damping ratio rj increases. 



By differentiating the both sides of Eq. ^ with respect to £, we get 

x(t) = — e~^(uj cos out — (3 sincjt), 
uj - ' 



x(t) 



v_o e -, t 

UJ 



(uj 2 — f3 2 ) sin uut + 2f3uj cos ut 



(8) 
(9) 



After solving Eq. ^ for t c in view of @, the following expression for the impact duration 
can be obtained [T3"] : 

1 -2Bu , . 

-Atan - (10) 



where the first positive value of the many- valued Atan function should be taken. 



Using properties of the Atan function, we rewrite Eq. (10) as follows [Uj : 

2/3uj 



tr = 



1 



id 



7T — atan 



co 2 -/3 2 
2f3u 

atan — — , u > p 



(3<u, 



(11) 



■ u 2 - f3 2 

Here, atan (z) is the principal branch of the arctangent function Atan (z). 



Finally, using properties of the atan function, we can rewrite formula (11) in a more 
simple form as 

2 u 

t c = — atan — . 12) 

u y ' 



Let 7] denote the loss factor, i. e., 



OJ 



(13) 
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Fig. 3. Viscoelastic Kelvin- Voigt impact model. Behavior of the main impact variables with time 
for the following values of the damping ratio: 77 = 0.1 (a), 77 = 0.3 (b), 77 = 0.5 (c). 



Then, Eq. (12) can be rewritten as 



at an 



(14) 



Recall that we assume that r] G [0, 1]. Also, note that in view of the notation 0, we have 

V=^^- (15) 



2\/km 



The velocity of the indenter at separation can be obtained by the substitution of (14) into 
^ in the following form: 



x(t c ) = -v exp(-/3i c ) 
= — v exp 



2r\ y/1 — rf 
atan 



Vi -n 



v 



(16) 
(17) 



From Eq. (16), it follows that the coefficient of restitution, e*, which is defined as the ratio 
of the velocity at separation |i(t c )| to the velocity of the indenter at incidence \x(0)\ = Vo, 
is given by 



f 2/3 uj 

e* = exp< atan — 

I uj p 



= exp < 



2r] y/1 — ri 2 
atan 



(18) 
(19) 



The peak value of the indenter penetration occurs at the instant t = t m , when x(t m ) = 0. 
In view of ([8]), we have 



t 



1 U) 

m — ~ atan — 
uj p 

1 



arcsin yl — r] 2 . 



(20) 
(21) 



Substituting the expression Q into Eq. @, we obtain the maximum penetration 
x(t m ) in the form 



^0 

Vo_ 

UJ 



exp 



exp 



f3 uo 

atan — 

uj p 

V 



arcsin J 1 — rf 



(22) 
(23) 



Note that from (14) and (21), it is readily seen that t m = t c /2. 



The peak value of the contact force F occurs at the instant t = t Mj when F(£m) — 0. 
According to Eqs. ([6]), ([8]), we obtain 



/■•(/) " ,r ° ' — ' 2 * 2 



UJ 



exp(— j3t) {uj — f3 ) sin ut + 2(3uj cos out 



f(l-2r] 2 ) I / \ 

= tyivqUJo exp(— rjuj t) — , sinu^y 1 — rft + 2r] cos uj^J 1 — rj 2 t . (24) 

V V 1 - T / 



Differentiating the previous expression, we can reduce the equation F(£m) — to the 
following one: 

1 - rj 2 (l - 4?7 2 ) cos^m^o^ 1 ~^ 2 ) + ~ 3 ) sin^M^oV 1 - 7 ? 2 ) = °- 



Thus, for 77 e (0, 0.5), we obtain 



_ 1 uj{uj 2 - 3f3 2 ) 
tM ~ u atan f3(3u 2 - W) 



1 



atan 



VT^p(l - 4ry 2 ) 
77(3 — 477 2 ) 



(25) 



For r] e (0.5, 1), the maximum value of the contact force F M = F(i M ) takes place at the 
initial instant t — 0. 



Substituting (25) into Eq. (24), we get 



V 



F M = mv u)o exp - 

F M = mv cu 2r], rj e (0.5, 1). 



atan 



71^2(1 — 4r7 2 )^ 
77(3 — 4?y 2 ) 



77 e (0,0.5), 



(26) 
(27) 



Note that the function F M {rf) defined by Eqs. (27) and (27) is continuously differentiate. 



Fig. [4^ shows the monotonic behavior of the dimensionless characteristic time moments 
tm/^0j t c / (2cl; ), ^m/^o with the damping ratio 77. Recall that t m = t c /2. The variations of 
the relative maximum contact force Fm/(™o^o) and displacement u)oX m /vo are presented 
in Fig. [4)3. It is interesting to observe the non monotonic behavior of Fm with the minimum 
at 77 « 0.26493. 



3 Viscoelastic Maxwell impact model 



Assuming that the cartilage layer's response to impact loading is modeled schematically 
as a serial combination of linear spring k and dashpot b (Fig. [5]). The force-displacement 
relation is given by the following differential equation [12]: 
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Fig. 4. Viscoelastic Kelvin- Voigt impact model. Behavior of the main impact parameters t m , t c , 
(a) and x m , Fm (b) with the damping ratio. 
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Fig. 5. Impact viscoelastic Maxwell model. 



From (28), it follows that 



F-*/«p(-J( t -r)}f(r)4r. 



(29) 



The differential equation of the impact mx+F = in view of (28) results in the third-order 



equation 



with the initial conditions 



... k k 

x + —x H x = 

b m 



x(0) = 0, x(0) = v , x(0) = 0. 



(30) 
(31) 



The impact problem (30), (31) has the following solution | [T3]fT5] : 

^ ^ 2Qv 



x(£) = — exp(— (cJot J \ sin uot — 1Q cos uot > + 



CJ 



x{t) = vq exp(— ((jj t) ^cos out + sincjt j. 



Here we used the notation 



k 



k 



(32) 
(33) 

(34) 



The variation of the contact force during the impact interaction is 

kv 



F = exp(— C^o^) sincji. 



(35) 



The impact duration t c is determined by the condition F =0. Thus, according to 



t=t c 



(35), the following relation takes place [T5IT5] : 



7T 



UJ 



(36) 



Substituting the value (36) into Eq. (33), one gets the coefficient of restitution in the form 

ttC 



e* = exp 



(37) 



Fig. [6] shows the behavior of the dimensionless quantities ojqx/vq, F/(m^ ^o), and x/vq 
with respect to time. Observe that the time moment £ m , when the indent er displacement 
reaches its maximum, approaches the final moment of impact as the damping ratio ( 



increases. 



According to Eq. (|33j), the peak value of the indenter penetration occurs at the instant 

(38) 



2uj\ 



1 H — arcsin ( . 



7T 



The substitution of the value (38) into Eqs. (32) and (35) gives the maximum penetration 
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Fig. 6. Viscoelastic Maxwell impact model. Behavior of the main impact variables with time for 
the following values of the damping ratio: ( = 0.1 (a), £ = 0.3 (b), ( = 0.5 (c). 



and the corresponding force 



F - =k -^ e M-2^( i+ l 3Ks ^» <4o) 



From Eq. (35), it follows that the peak value Fm of the contact force occurs at the instant 

i 



tu — — atan 



c 



(41) 



Substituting (41) into Eqs. (35) and (32), we obtain the maximum contact force 



„ kv C 
F M = exp<^ 



atan 



c 



and the corresponding displacement 



atan 



(42) 



(43) 
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Fig. 7. Viscoelastic Maxwell impact model. Behavior of the main impact parameters t mi t c , tM 
(a) and x m , Fm (b) with the damping ratio. 

Fig. [7] shows the monotonic behavior of the dimensionless characteristic time moments 
tm/uo, t c j (2cjo), ^m/^o with the damping ratio (. The variations of the relative maximum 
contact force F M / '(mvoUJo) and displacement ujQX m /vQ are presented in Fig. [7)3. 



Finally, as it was observed [13] , although certain quantities of the Maxwell impact model 
are equivalent to the so-called half-period Kelvin -Voigt impact model, the inherent 
physics of these models are completely different. 



4 Standard solid model 



There are two schematic representations of the standard linear solid model (Figs. [8] and 
9b. The force-displacement relationship is given by the following two equations: 



(fci + k 2 )F + bF = k\k 2 x + k\bx, 
H\F + f3F xix 2 x + $(k\ + x 2 )x. 



(44) 
(45) 



K 



Fig. 8. Standard solid model. Configuration based on the Kelvin- Voigt model. 
The instantaneous and long-term moduli are 

hk 2 



ko = fa = xi + x 2 , = xi 



fa + k 2 ' 



(46) 
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Fig. 9. Standard solid model. Configuration based on the Maxwell model. 
The relaxation time is equal to 

b 



fa + k 2 x 2 ' 



(47) 



The differential equations (44) and (45) are equivalent to the force-displacement relation- 
ship 

dx 



F= [ k(t-r)^(r)dr 
J dr 



(48) 



with the relaxation stiffness 



k(t) = k QO + (k - fcoo) exp 



Let us also introduce the notation 



k 
k 



Note that p e (0, 1) in view of ([46]). 

The differential equation of impact 

mi + F = 0, te [0,* c ], 

where the contact force F is determined by Eq. ([44]) , can be written as 

... , (h + k 2 ) k x k x k 2 

x + -x H x H —x 0. 

b m mb 

By introducing the non-dimensional time 



tr 



Eq. (52) can be reduced to the following equation: 

x m + x" + Ax' + Apx = 0. 



Here we introduced the notation 



A = 



m 



The initial conditions for Eq. (|52j) are as follows: 

x(0) = 0, x(0) = v , x(0) = 0. 



The solution to the problem ( 52 ) , ( 56 ) is given by the following formula 



x(t) 



CiPi-A^ + d 2 ] 



[(1-ftXAx + sin 



, ^ , x (it) ( Pit 

— Ai) cos — > exp 



1 - Ai)r^ 



r R ; 1 (a - Ai) 2 + Ci exp 



Ait 



(49) 



(50) 



(51) 



(52) 



(53) 

(54) 
(55) 

(56) 



(57) 



Here, — Ai and —(/3i ± are the roots of the characteristic equation corresponding to 

(58) 



Eq. (54). In other words, the following factorization takes place: 

z 3 + z 2 + Az + Ap = (z + X^iz 2 + 2(3 lZ + ft 2 + Ci). 



The discriminant of the characteristic equation is 

D = 4A(A 2 + p) - A 2 (l + 18p - 27p 2 ). 



(59) 



We underline that formula (57) holds true when D > 0. In this case, we have 



Ai 

Ci 

where 



3 3 


(1 - 3A) 


3d ' 


1 Ci 


(1 - 3A) 


3 6 


6Ci ' 




- 3A) 


6 


6Ci 



Ci = OQi + 2-9A + 27Ap), Qi 



(2 - 9A + 27Ap) 2 - 4(1 - 3Af 



5 Perturbation of the Kelvin— Voigt model 



Taking into account (46) and (50), we rewrite Eq. (44) in the following form: 

kocF + p(l - p)bF = k^x + (1 - p)k OQ bx. 



(60) 



Now, letting p 0, we arrive at the equation 

F kr^rX + 6x, 



(61) 



which coincides with Eq. ([!]). Thus, for small values of p, the standard solid model (60) 
is a perturbation of the Kelvin- Voigt model (61). 



Let us introduce the notation 



2 ^oo Q b ft 



m 



2m 



UJ 



(62) 



Then, the parameters (47) and (55) can be evaluated as 



r R = 2 V p(l - p) — , A = 477 2 p(l " P) 2 - 



(63) 



In view of (63), the discriminant (59) and the roots of the characteristic equation (58) 



can be asymptotically evaluated as follows: 



D = 16i] 2 (l -ri 2 )p 2 + 0(p 3 ), p^0, 
A 1 = l -4?7 2 p + 0(p 2 ), 

A = 2?y 2 p + 0(p 2 ), Ci = 2 V Jl - rj 2 P + 0(p 2 ). 




Fig. 10. Perturbation of the Kelvin- Voigt model. Relative errors of the asymptotic formulas (64) 



and (65) for the duration of impact (a) and the coefficient of restitution (b). 



The accuracy of the asymptotic approximations (64) and (65) is presented in Fig. 10 



Note that the asymptotic formulas (64) and (65) are not uniformly valid as rj 1 



6 Perturbation of the Maxwell model 



Now, taking into account (46) and ( |50| ), we rewrite Eq. ( |44| ) as follows: 

k F + (1 - p)bF = pklx + (1 - p)k bx. 



(66) 



Again, by letting p — >■ 0, we obtain the limit equation 

F F 



k 



(67) 



which coincides with Eq. (28). Thus, for small values of p, the standard solid model (66) 



can be regarded as a perturbation of the Maxwell model (67) 
Let us introduce the notation 



k 



m 



c = 



2uj b 



(68) 



In view of (68), the parameters (47) and (55) can be evaluated as 

(i - pf 



TR 



1-p 



A 



4C 2 



(69) 



characteristic equation (58) as follows: 



Now, taking into account (69), we expand the discriminant (59) and the roots of the 



Ci 



16C 6 8C 6 
+ 0(ft 2 ), fr = \- 

vr^e vw^(2 + c 2 ) 



A 1= p + 0(p 2 ), Pi = \-\ + 0{p\ 



2C 



- p- 



2C(1 " C 2 ) 



+ 0(p 2 ). 



Consequently, we obtain the following asymptotic approximations for the impact duration, 
i c , and the coefficient of restitution, e*: 



u)ot c 



7i 2vrpC 



exp 



ttC 



+w 1+ 1 



v'R (i-C 2 ) 3/2: 



2(1 — C 2 ) 3/2 , 



exp 



ttC 



(70) 
(71) 



The accuracy of the asymptotic approximations (70) and (71) is presented in Fig. 11 



Note that the asymptotic formulas (70) and (71) are not uniformly valid as rj 1. 



7 Drop weight impact. Viscoelastic Kelvin— Voigt model 



Due to Newton's second law, the differential equation of the drop weight impact has the 
form 

mx + bx + kx = mg ) t G [0, t c ] , (72) 
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Fig. 11. Perturbation of the Maxwell model. Relative errors of the asymptotic formulas (70) and 



(71) for the duration of impact (a) and the coefficient of restitution (b). 



where g is the gravitational acceleration. 



The initial conditions for Eq. (72) are 

x(0) = 0, x(0) = v . 



(73) 



The drop weight impact problem (72), (73) has the following solution: 



x 



(t) = JL(\- e-^coscjt) + -(v - ^]e- 0t smcut, 
tot v ' oj \ ujA J 



(74) 



x(t) = v e~ pt coscjt + — — ^-e~ 0t smut. 



(75) 



Here we used the notation ([7]). 



According to Eqs. (74), (|75|), the reaction force F(x, x) = kx + bx is given by 



= 6:q < 1 + e p , sin uot — cos cji 



H^^+tBw™^)' (76) 



where we introduced the notation 



eo = —. (77) 

UJqVq 



The problem ( 72 ) , ( 73 ) was studied in [16] , where the existence of the parameter domain 
of "plastic impact" was established. This means that for any r] > 0, there exists a unique 
value of £q such that for all e > £q we have x) > in the time interval t G (0, +00). 
The critical value £q of the parameter £ determines the critical value Vq of the initial 
velocity v below which there is no rebound effect. 

With the aim of application to the drop weight impact testing, we consider the problem 



(72), (73) for small values of the dimensionless parameter £ G and construct an asymptotic 



solution for the coefficient of restitution. 



Let and e° be the impact duration and the coefficient of restitution for the Kelvin- Voigt 



impact model, correspondingly. According to Eqs. (12) and (19), we have 



t° = 



atan 



e* = exp< 



2rj 



VI -V' 



atan 



(78) 



Now, solving the transcendental equation x] 
terms of the first order inclusive, we obtain 



t=t c 



by a perturbation method to 



(79) 



e* ^ e°(l - 2e r]), 



(80) 



where t® and are given by Eqs. (78). 



From the asymptotic formulas (79) and (80), it is clearly seen that the gravitational effect 



increases the duration of the impact process and decreases the coefficient of restitution. 
But it is more interesting to observe that the coefficient of restitution e* increases with 
velocity v 0) since the parameter £0 is inversely proportional to v . That is why the effect 
of decrease in the coefficient of restitution in the drop weight impact test experimentally 
observed in [7] for the velocity range vq G (0.7,1.4) m/s and extrapolated for the low 
velocity region by means of the nonlinear Kelvin- Voigt model F(x : x) = kx + c\x\x with 
no account for the impactor weight cannot be explained by the linear viscoelastic Kelvin- 
Voigt model considered in this section. 



8 Drop weight impact. Viscoelastic Maxwell model 



By applying the approach [15J, the differential equation of motion mx + F = mg with the 
initial conditions x(0) = and x(0) = v in view of the constitutive relationship (28) can 
be reduced to the following problem: 

k .. k . ko ^ > v 

x + -x-\ x--=- = 0, (81) 

b m b 

x(0) = 0, x(0) = v , x (0) = p. (82) 
The drop weight impact problem plj), ((821 has the following exact solution: 



x(t) 



— e'^A — (1 - 2C 2 - £ C) sin wt - (2( + e ) coswf. 



(83) 



±(t) = v oe - Cwo *|^(C + £ ) sinwt + coswr. 
+ 2C£ ^o- 



(84) 



Here we used the notation (34), (77). 



According to Eqs. (83) and (84), the reaction force F = mg — x is given by 



mvouuo 



= — e ^ ot \ (1 + £ C) sincjt — —£o cos out > + e 

U) { U)Q J 



(85) 



Now, let t® and e° be the impact duration and the coefficient of restitution for the Maxwell 



impact model, correspondingly. According to Eqs. (36) and (37), we have 



7T 



exp 



(86) 



Applying a perturbation method, we find 



2(e , 



(87) 



where t£ and e® are given by Eqs. (86). 



It is interesting to note that the asymptotic formulas (79) and (87) coincide. Clearly, 



the same conclusions can be drawn about the influence of the gravitational effect in 
the framework of the viscoelastic Maxwell drop weight impact model as those that were 
formulated in Section [3 



9 Short-time asymptotic solution of the indentation problem for a thin tripha- 
sic layer 



We assume that the deformational behavior of articular cartilage is modeled in the frame- 
work of linear biphasic theory [ID], which represents the biological tissue as a mixture 
consisting of a porous solid phase and a fluid phase (mobile interstitial water). The con- 
stitution equations for the solid and fluid phase stresses, cr s and cr^ are given by 

a s = -^ s pl + A s tr(s)I + 2//s, a f = -cj) f pl. 

Here, $ is the fluid volume fraction (porosity), (j) s = 1 — $ is the solid volume fraction, p 
is the true pressure of the fluid, X s and ji s are the Lame constants, which together define 
the aggregate modulus Ha = A s + 2/x s , e is the strain tensor of the solid phase, and I is 
the identity tensor. Note that the fluid phase is assumed to be intrinsically incompressible 
and inviscid. 

The continuity equation for the mixture and the momentum equations for each phase are 
given by 

div(0 s v s + <f> f v f ) = 0, 

div<r s + ^^(v' - v s ) = 0, diva f - ^^(v f - v s ) = 0, 

where v s and are the solid and fluid velocities, respectively, and k is the permeability 
of the solid phase. 

Let us consider an axisymmetric contact problem for a thin biphasic layer indented with- 
out friction by a rigid impermeable cylindrical indenter. It is assumed that the contact 
radius a is much larger than the cartilage layer thickness h (i.e., h/a <C 1). In this case, 
according to [17] . the vertical displacements w(r,t) of the boundary points of the articu- 
lar cartilage tissue at the contact zone can be approximated by the following asymptotic 
formula: 

, N h 3 { 1 d ( dP, A fc^ [19/ dP , A 7 ] 







Here, P(r,t) is the contact pressure. It is assumed that the cartilage layer is bonded to a 
rigid impermeable substrate, that is there is no solid displacement at the cartilage-bone 
interface and no fluid flow through the bone [17] . 



In view of (89), the contact condition that the boundary points of the cartilage layer 
acquire a constant vertical displacement —So(t) (due to the action of the indenter) can be 
written as 

w(r,t) = -6 (t), r<a. (90) 



The substitution of (89) into Eq. (90) results in an integro-differential equation 



lfl ^))^/;l(^))*-^ < 91 > 



r dr\ dr 



which requires imposing a suitable boundary condition at the edge of the contact zone, 
i.e., at r = a . 



In order to impose the mentioned boundary condition, we note that at the initial moment 
of contact t — 0, formula (89) simplifies as follows: 



w(r, 0) 



h 3 i d 



3/i s 3r dr V dr 



dP 

r^-—(r, 0) 



(92) 



Comparing formula (92) with the known asymptotic solutions for thin elastic layers 



[T8lll9120j . we conclude that the instantaneous deformational response of a thin bipha- 
sic layer coincides with the response of a thin bonded incompressible elastic layer. Thus, 
by this analogy, we will require that P(r, t) —± as r — » a, that is the contact pressure is 
assumed to vanish at the edge of the contact area. 



As a result of integration of Eq. (91) with respect to the radial coordinate, we arrive at 
the following integral equation: 



P(r, t) + 3 -jff P(r, r) dr = |j W(a 2 - r 2 ). 



(93) 



Now, in order to derive the relationship between the indenter displacement and the contact 
force 



F(t) = 2tt J P{r,t)rdr, 



we multiply both sides of Eq. (93) by 2nr and after that we integrate the equation obtained 
with respect to r from to a. As a results of this operation, we get 



F(t) 



3fl s K 



(94) 



Further, by inverting the Volterra integral operator on the right-hand side of Eq. (94), we 
obtain 



(95) 



where we introduced the shorthand notation 



X 



3/I s K 



(96) 



Finally, after integrating by parts, Eq. (95) yields 



F(t) 



3/U s a 4 
16/i 3 



)(0) + / - 



-X(*-r)^( T)dT \ 

dr 



(97) 



In impact problems, under the assumption that 

So(0) = 0, 

the force-displacement relationship (l97|) takes the form 



(98) 



Here we introduced the notation tr = 1/x- In view of (96), we have 

h 2 

TR 



(99) 



while comparing (99) with (29), we get the stiffness coefficient 

3/x s a 4 



k 



16h 3 ' 



(100) 



It should be emphasized that Eq. (98) represents a short-time asymptotic approximation, 
which is valid for moments of time t such that HA^t/h 2 <C 1. For typical human cartilage 
material properties, Ha = 0.5 MPa and k = 2 x 10 -15 m 4 /Ns. Thus, assuming a typical 



cartilage thickness h — 1 mm, we get h 2 /(Hak) = 10 3 s; thus, the asymptotic model (98) 



certainly remains valid for up to 100 s, which is well in the range of usual values of impact 
durations. 



Comparing Eq. (98) with Eq. (29), we see that the short-time deformational response of 



a thin biphasic layer bonded to a rigid impermeable substrate under the action of a fric- 
tionless flat-ended indenter is mathematically equivalent to that of a thin incompressible 
layer following the Maxwell viscoelastic model. Note that the Maxwell's model based per- 
turbation model considered in Section [6] could be useful for modeling the impact response 
of articular cartilage (or artificial tissues for its replacement) in the whole time range, i.e. 
in the short-, medium- and long-time range. 

We also emphasize that the biphasic model is not equivalent to a viscoelastic model, 
because the biomechanical response of a poroelastic material such as articular cartilage is 
crucially dependent on the boundary conditions for the sample. In particular, viscoelastic 
equivalents of the deformational response of an articular cartilage sample subjected to the 
same simple loading protocols in confined and unconfined conditions will be essentially 
different, especially in the short-time range. Thus, in comparing experimental results from 
different sources, a particular attention should be paid to the fixation conditions for tissue 
samples. 



10 Key features of non-linear impact 



To illustrate the application of the developed linear theory of viscoelastic impact, let us 
analyze the experimental data obtained in [5j for drop-weight impact testing (with the 



impactor mass m = 100 g) of isolated bovine articular cartilage samples of 5 mm diameter 
(correspondingly, the radius of the samples is a = 2.5 mm). In [5], the force data, F(t), 
were converted to engineering stress, cr(i), by dividing them by the original cross-section 
area of the sample, 7ra 2 , i. e., 

a W = — 2' 



The effective strain, e(t), was evaluated by dividing the measured impactor displacement, 
x{t) ) by the sample thickness, which is assumed to be 0.5 ± 0.11 mm, as follows: 



e(t) 



x(t) 



(Here, stress and strain are assumed to be positive in compression.) The stress-strain 
relationship was differentiated to obtain the incremental dynamic modulus 



E, 



dyn 



da 
de 



The maximum incremental dynamic modulus, i? ma x, was found, and the modulus Ei at 
stresses of 10 MPa was determined to enable comparison of dynamic moduli at constant 
value of stress. The initial impact velocity was calculated from the drop height, /i , by the 
well-known formula vq = yj2gh$. 



The incremental dynamic modulus can be evaluated as a function of time in the form 

(101) 



p (A h F(t) 



e(t) na 2 x(t) 

In the case of the Maxwell model (see, Section [3]), we will have 

F(t) ^cos out — ((ujq/uj) sin ujt 
x(t) cos out + ((ujq/uj) smut ' 

where k is the stiffness coefficient. 



(102) 



First of all, observe that in view of (101) and (102), the variation of Ed yn (t) does not 
depend on the impact velocity v . In other words, the time variation of the incremental 
dynamic stiffness in the linear viscoelastic impact tests remains the same for different 
initial impact velocities. We emphasize that this conclusion is valid for a general linear 
viscoelastic law. Second, from (101) and (102), it follows that the value of E^ yn (t) gradually 
decreases to zero with increasing contact force F(t) (when F(t) > 0). Thus, we arrive at 
the formula 

£max = £dyn(0). (103) 



Further, in order to evaluate Eio, we need first solve the equation 



F(t 



10 



ira oio, 



where crio = 10 MPa. In view of (35), Eq. (104) takes the form 



exp(- C^o*io) sin ut 10 = 



7TCL OJCTio 



kv { 







(104) 



(105) 



Here, £, cjq, and are independent of vq, and are determined by formulas (34). 



Now, from (105), it is seen that the value of the time moment tio depends on the initial 
velocity v . Thus, the Maxwell impact model (and generally speaking, any linear vis- 
coelastic model of impact) predicts that the value of Ei increases with increasing impact 
velocity v . 

Table 1 

Impact parameters for isolated bovine articular cartilage samples [5] 
h (mm) v (m/s) £ max (MPa) E 10 (MPa) cr max (MPa) e max e* Am (%) 



25 


0.70 


86 ±22 


75 ± 13 


15.6 ±2.9 


0.48 ± 0.06 


0.64 ±0.08 


2.2 


50 


0.99 


100 ± 32 


71 ± 16 


24.5 ±3.5 


0.60 ±0.13 


0.46 ±0.14 


2.5 


80 


1.25 


118 ±33 


73 ± 12 


34.2 ±5.0 


0.62 ±0.11 


0.47 ±0.05 


5.7 


100 


1.40 


128 ± 28 


72 ± 13 


40.5 ±4.6 


0.68 ± 0.09 


0.41 ±0.08 


9.9 



Table [T] shows that the impact testing [5] was performed in the non-linear regime with 
maximum compressive strains of 50-60%. That is why, the prediction of the linear impact 
model concerning i? max are not fulfilled. Furthermore, the linear theories of impact predict 
that the maximum contact force F M (correspondingly, the maximum contact stress cr max = 
Fm/^o 2 )) and the maximum displacement x m (correspondingly, the maximum strain 
e max = x m /h) are proportional to vq. On the other hand, the data from Table 1 show 
that the ratio cr max /e max increases with increasing v . This fact also clearly indicates the 
non-linear deformational behavior of cartilage at high level of strain. Note here that the 
ratio cr max /e max is raited to the so-called pulsatile dynamic modulus (see, in particular, 

EH). 

Concerning the coefficient of restitution e* note that it is not constant, as it would be if 
the cartilage deformation were described by the Maxwell model (see formula (|37|)). 



Finally, the last column of Table [T] gives the values of percentage increase in mass of 
each sample after 24 h immersed in PBS following impact loading. This is indicative of 
increasing amounts of damage in the cartilage samples [5] . 



11 Discussion 



Consider now the general case of linear viscoelastic force-displacement relationship 

F = [ k(t-s)^(s)ds (106) 
J as 
o 

with the relaxation stiffness 

Here, ko is the initial stiffness, tr is the characteristic relaxation time, \P(t) is the dimen- 
sionless relaxation function with r being a dimensionless independent time-like variable. 



Making use of the change of variables 

t = t r t, x = v t r £, 



(107) 



we transform the impact equation mx + F = and the initial conditions x(0) = 0, x = v 
into the following problem: 



f + aj V(t-<t)^{<t) da = 0, 



(108) 



£(0) = 0, £'(0) = 1. (109) 
Here prime denotes differentiation with respect to t, and we introduced the notation 



a = 



m 



(110) 



Note that for the Maxwell model (see Section [3j Eq. (29)) we have k = k, tr = b/k, and 
« = 1/(4C 2 )- 

Furthermore, according to Eqs. ( 107| )), the variable impact velocity is 

±(t)=V £'(T). 



Let r c be the dimensionless duration of the impact process. Then, the coefficient of resti- 
tution can be found as 

e„ = -e'(r c ). (Ill) 

From Eqs. ( 108| )) and ( 109| )), it is evident that r c is a function of a only and does not 
depend on v Q . Thus, in view of ( HID ), we conclude that the coefficient of restitution e* is 
constant with respect to the initial impact velocity vq. 

It can be shown that the same qualitative conclusions are drawn from the linear biphasic 
model [10J for articular cartilage deformation. In this case, the parameter t#, which enters 
Eqs. (107)), can be defined as a typical diffusion time = h 2 /{K J HA) ) where h is the 
cartilage layer thickness, k is the cartilage permeability, and Ha is the aggregate modulus. 

Remark 1 Let us consider the question of applicability of the coefficient of restitution 
for diagnosis of the state of health of the tissue. In the framework of the asymptotic model 
(98), according to Eq. (31), we will have 



e* exp 



(112) 



where (see the last formula (34)) 



c 



k 



2u)ob 



with k and b defined as follows (see Eq. (99) and (100\) ): 

3/x 5 a 4 b h 2 
16/i 3 ' k R 3/jl s k 



(113) 



(114) 



From (112) it is readily seen that the coefficient of restitution e* decreases with increasing 
loss factor (. At the same time, in view of (113) and (114), we have 



It is known [22] that for articular cartilage in the early stages of osteoarthritis, the fol- 
lowing degenerative changes are observed: increased permeability, k, increased thickness 
of the cartilage layer, h, reduced shear modulus, ji s , and/ or a combination of these ef- 



fects. Formula (115) implies that increasing the permeability of the cartilage results in a 
decrease of e*, while increasing the cartilage thickness and decreasing the shear modulus 
(both create a softening effect) apparently results in an increase of e*. Because it is known 
JEEjj that osteoarthritic cartilage may show a dramatic (more than 6-fold) increase in the 
hydraulic permeability k, it can be expected that the overall change in the coefficient of 
restitution e* will be negative. 



It should be also observed that formulas (115) and (31) imply an increase of the coefficient 



of restitution e* with increasing the cartilage thickness h, whereas the experimental data 
presented in Fig. [7] apparently show an inverse tendency. 

Remark 2 It is known [4,5] that impact loading of articular cartilage at high impact 
stresses typically result in fissuring of the articular cartilage surface. At the same time, 
the formation of cracks allows to absorb greater amounts of energy as well as dramatically 
affect the deformation resistance of cartilage resulting in change of the parameters of the 
impact model. In other words, the mechanical properties of the tissue do not remain the 
same to the end of the impact process. 

Observe that the biphasic theory incorporating Lame parameters assumes that the ma- 
terial of solid phase is linearly elastic in order for these to have unique values. But if the 
material is viscoelastic these parameters are difficult to define and they become functions 
of deformation and/or time, if they are meaningful at all. Furthermore, there is an intrinsic 
circularity problem associated with using the aggregate modulus Ha, which is evaluated 
at equilibrium after the interstitial water is squeezed out, to define the mechanical prop- 
erties that are then assumed to pertain during the impact deformation. Thus, the fact 
that the biphasic theory provides a good fit to measured curves in the creep and stress 
relaxation tests can be basically considered as a consequence of a curve-fitting procedure 
with a minimum of three free parameters rather than a derivation from first principles. 
In other words, it remains to be an open question on the efficiency of mixture models for 
articular cartilage at high strain rates. 

In the present study we addressed the question of whether the main features of articular 
impact observed in [6117] could be qualitatively predicted using a linear viscoelastic theory 
or the linear biphasic theory. It is to note that the deformations encountered in impact 
tests should be small enough for the linear theories to apply. With respect to engineering 
polymers note that the linear theory of viscoelasticity may hold reasonably well even up 
to some 5 — 10% extension, in particular for certain rubbers [24J. 



12 Conclusions 



The results of this study based on the linear viscoelasticity imply the following properties 
of the linear impact models: 

1. The coefficient of restitution e* is a function of the damping ratio ( alone. This means 
that e* does not depend on the impact velocity Vq, but it depends on the impactor mass 
m and the sample thickness (through the stiffness k). 

2. The impact duration t c is inversely proportional to cjq, that is t c is proportional to \fra^ 
and depends on the damping ratio as well. The impact duration does not depend on the 
impact velocity v . 

3. The maximum displacement, x m) and the maximum contact force, F M) are proportional 
to ^Jm and v . 

4. The time variation of the incremental dynamic stiffness dF/dx remains the same for 
different initial impact velocities. 

5. In the drop weight impact test, the gravitational effect increases the impact duration 
t c and decreases the coefficient of restitution e*. At that, the coefficient of restitution 
increases with the impact velocity v . 
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